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We discuss an equivalence between the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formula- 
tion of the Einstein evolution equations, a subfamily of the Kidder-Scheel-Teukolsky formulation, 
and other strongly or symmetric hyperbolic first order systems with fixed shift and densitized lapse. 
This allows us to show under which conditions the BSSN system is, in a sense to be discussed, hy- 
perbolic. This desirable property may account in part for the empirically observed better behavior 
of the BSSN formulation in numerical evolutions involving black holes. 
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I. INTRODUCTION 

o : 

When numerically integrating the Einstein equations, it is common to split them into constraint and evolution 
equations. Most attempts usually solve the "free evolution" problem, solving the constraints only at the level of 
^"v initial data. As in all gauge theories, the evolution equations of general relativity are not uniquely defined. One can 
redefine them by adding constraints, since the latter vanish on solutions of the whole system of equations. It has been 
increasingly recognized in the last few years (see jy for a list of references) , that it is desirable to take advantage of 
this freedom to make the evolution equations strongly- (and preferably symmetric-) hyperbolic. There is a body of 
mathematical literature ^ that one can draw from to prove stability (in the sense of Lax's theorem) for discretization 
schemes based on well posed formulations. In hyperbolic formulations one also has better knowledge of the modes of 
, the equations, something that is of interest at the time of specifying boundary conditions that are consistent with the 
■ constraints and the evolution equations || [I| . 

The most traditional of the evolution schemes, the Arnowitt-Deser-Misner (ADM) scheme, when written in first 
order form, is only weakly hyperbolic (WH) Q. There has always been a feeling that this accounted in part for the 
difficulties encountered trying to evolve black holes and other spacetimes with this evolution scheme. More recently 
||, it was rigorously shown that popular discretizations of the ADM equations in particular, and WH formulations 
in general, simply do not lead to convergent schemes. 

It was first recognized by Shibata and Nakamura Q| and further elaborated by Baumgarte and Shapiro Q that 
if one performed a conformal-traceless (CT) decomposition of the ADM variables, the resulting scheme (commonly 
known as BSSN), at least in certain situations, appeared to have less numerical difficulties than the ADM system. 
Also, anecdotal evidence suggested that slightly different implementations of the BSSN scheme seemed to perform 
^Jy quite differently in numerical simulations. 

In this paper we would like to point out that the BSSN scheme, under certain conditions, is a reduction (from first 
to second order in space) of some strongly or symmetric hyperbolic systems. This may account in part for its better 
, performance in numerical simulations. The conditions under which these reductions can be established depends rather 
1 subtly on details of the implementation of the BSSN scheme, this may account for why some implementations appear 
| better behaved than others. 

We will show that the BSSN scheme and its generalization by Alcubierre et. al. [Q can be mapped to the Kidder- 
Scheel-Teukolsky (KST) Q system of evolution equations in the linearized regime. The latter is a multi-parameter 
family of evolution schemes, and we will discuss for which values of the parameters the correspondence holds. The 
BSSN system is a first order in time and second order in space formulation, whereas the KST system is first order in 
time and space. Therefore we will need to perform an enlargement of the BSSN system to establish the correspondence. 
We will discuss in detail the precise nature of the correspondence below. We will show that a similar enlargement 
exists at the nonlinear level, which does not yield the KST system but which does yield a strongly or symmetric 
hyperbolic one. The latter system is closely related to one introduced by Friedrich and Rendall [jio] . 

We will start by presenting, in Section II, the KST system and summarizing their main properties. Next we present 
a CT decomposition of this system that shares its level of hyperbolicity. We then discuss in detail the correspondence 
in linearized gravity between this CT decomposition of the KST system and the BSSN equations. In Section III we 
show that by introducing extra variables to make the BSSN system first order in time and space one can get strongly 
and symmetric hyperbolic systems provided that the lapse is appropriately densitized. We end with a discussion of 
the relevance of these results for numerical simulations. 
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II. THE BSSN-KST CORRESPONDENCE 



In this section, we show that for the case of linearization around flat space-time in Cartesian coordinates the BSSN 
equations are equivalent (in a sense to be defined below) to the linearized KST equations with an appropriate choice 
of parameters. As a byproduct of our construction, we present a CT decomposition of the KST system that shares 
its same level of hyperbolicity, even at the nonlinear level. 



A. The KST family of equations 



Starting from the ADM equations, KST derive a family of strongly and symmetric hyperbolic first-order evolution 
equations for the three- metric Qfo), the extrinsic curvature (-fTy), and the spatial derivatives of the three-metric 
(dkij = dhQij), which generalizes previous well posed formulations ]Tl| , jL2| . A priori prescribing the densitized lapse 
exp (Q) as a function of spacetime, the lapse N is given by 

N = g°e { >, (I) 

where g is the determinant of the three-metric and a a parameter. The shift vector, is also assumed to be prescribed 
as a function of spacetime. Here we use the same notation as in KST except that we define the derivative operator 
along the normal to the hypersurfaces as 

do = jj(dt - £p). 
The vacuum evolution equations have the form 

dogij = -2Kij , (2) 
d K i:i = ~d k d kij + - 2 if <l + ig ab %d|afcb) ~ Q + CT ) 9 ai d {i d 3)ab 

+ Cg ab C a{lj)b +i gi] C - e-Qdid^e®) + TUj , (3) 
dodkij = -2dkKij + , qg k ( i C j ) + XQijCk + T^kij , (4) 

where the constraint variables 

C — \g ab d k {d a bk — dkab) + 7t, (Hamiltonian constraint) 

Cj = d a K a j — g ab djK a b + TZj, (momentum constraint) 

Ckij = d kij - dk9ij , (definition of d kij ) 

Cikij = d[id k ]ij ■ ("closedness" of d kij ) 

have been added to the right-hand side of the evolution equations with some free parameters £, 7, 77, \- Here, d k = g kl di 
and bj , d k are defined as the traces of dkij ■ 

bj = g kl dkij , dk — g^dkij ■ 

The "remainder" terms 1Z, IZi, IZij and IZkij are homogeneous polynomials of degree 2 in dkij, diQ and Kij, where 
the coefficients may depend on gij . Finally, the Lie derivative of the symbols dkij is 

£pd k ij = P l did kij + du 3 dkP l + 2d k i( l dj ) 0- + 2gi (l d 3) d k 0- ■ 

The evolution system has characteristic speeds {0, ±1, zt-^/Ai, ±\/A2, dzy/X^}, where 

Ai = 2(7, 

A 2 = 5(1 + 0*7 + 7(2 -»7 + 2x), 
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It should be noted that these are the characteristic speeds with respect to the do operator. This means that A = 0, 1 
correspond to propagation along the normal to the hypersurfaces or the light cone, respectively (physical speeds). 
The characteristic speeds with respect to the dt operator are obtained from these after the transformation 

(ih Nfi + f3 l n % , 

where n l is the direction of the corresponding characteristic mode. 

The conditions under which the system is completely ill posed (CIP), WH or strongly hyperbolic (SH) were found 
by KST. The system is CIP if any of the above speeds is complex, while it is WH if Xj > for j = 1,2,3 but one of 
the conditions (|5|, ^|) below is violated. For example, if the parameters (£, 7, rj, \) are zero the dynamics is equivalent 
to the ADM equations written in first order form with fixed densitizcd lapse and fixed shift (which are WH) . If a = 
as well, then the system is equivalent to the ADM equations with fixed lapse and shift (which are also WH). 

Finally, the system is SH provided that 

A;>0, for j = 1,2, 3, (5) 

A 3 = ^(3A! + 1) ifA 1= A 2 . (6) 

Here by SH it is meant that the principal part has a complete set of eigenvectors with real eigenvalues. Provided 
there is a smooth and uniformly bounded symmetrizer, this ensures well posedness Qj. 

In KST seven extra parameters are introduced and used to make changes of variables in Kij and dkij ■ When 
performing this change of variables the constraint Ckij — is also used in order to trade spatial derivatives of the 
three-metric for dkij ■ Thus, the equations with the new variables have different solutions off the constraint surface. 
However, one can see that at the linear level this change of variables does not involve the addition of constraints. 



B. Conformal trace-less decomposition of the KST system 



We start by performing a CT of the KST variables: 

9ij — 6 ^9ij ? 



K, 



ikij 



^ (A tj + ^k) 



2e klj + -cjk(i^j) — -ffyTfe + -gijdk 



Here has unit determinant, Aij is trace-less and e k ij is trace-less over all pairs of indices, and Tj is defined by 
Tj = bj — dj/3. In terms of the Christoffel symbols T\ 3 - belonging to the conformal metric cjij, Tj can also be expressed 
as (d = g lk dk ) 

Tj = g kl Tjki = d % gij , 

which is the additional variable used in the BSSN system (actually, in the original version by Shibata and Nakamura, 
the variable Fj — S kl dkgij is used instead of Tj, in the version by Baumgarte and Shapiro, r fc = g ^Tj is used.) 

With respect to this decomposition, the KST system (H-f|) becomes, up to lower order terms ("l.o.") quadratic in 
dkij, diQ and Kij, 

1 K, (7) 



<9o0 
9o9ij 
d K 

dod k 
doTj 

doekij 



G 

-2Ai 

1 + ^7 ) d k T k - ( ~- + 7 + o ) d k d k 



- Q d k d k (e Q )+l.o. 



-d k e kll + (1 + C)d k e {lj)k + e-^ ^(5 - 9()d (l T 3) - Q + a^j d^d j} - e-^(e«) 



TF 



(r? + 3 X )d s A sk --( v + 3 X + 3) d k K + (l.o.) k , 

ft-*) d s A S3 -^ v d 3 K+(l.o.) 3 , 

3 1 

-d k Aij + -g k ( t d s A j)s - -g l0 d s A sk + {l.o.) kl3 



(8) 
(9) 

+ (lo.)ij , (10) 

(11) 
(12) 

(13) 
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where the superscript TF denotes the trace-free part. Here, for simplicity, we have also assumed that the shift 
vanishes. However, a generalization in the presence of shift can also be worked out. 

Notice that in the derivation of Eqs. (|7 - 13 ) we have used the constraints Ckij = in order to re-express first order 
spatial derivatives of <?.y in terms of dkij in a way that ensures that the principal part of the transformed equations 
has the same spectrum as the original KST system. Therefore, for a given choice of parameters both systems have 
the same level of hyperbolicity. However, in general the dynamics of the evolution Eqs. ([7[jl3|) is not identical to the 
dynamics of the KST Eqs. (Q-Q) if we are off the constraint surface Ckij = 0. A similar situation occurs in KST when 
the variables Kij and dkij are redefined through a seven-parameter family of change of variables. 



C. Correspondence 



We start by looking at the linearized equations around flat spacetime, written in Cartesian coordinates. One can see 
that in this case all lower order terms ("Lo.") vanish, and the two systems (||^|) and (0-|l^) are dynamically equivalent 
(i.e. have the same set of solutions, even off the constraint surface). In the following, we will call the linearized version 
of equations ([^1^) system A. 

Since the BSSN system has second order spatial derivatives and the KST system has only first order spatial 
derivatives, care needs to be exercised concerning the constraints that materialize the correspondence of the spatial 
derivatives. Basically, we need to show that if we start with initial data that is on the submanifold of constraints that 
relate the second spatial derivatives to first spatial derivatives of extra variables (all solutions of the BSSN system 
automatically lie on this submanifold), we automatically remain on it. Otherwise, no correspondence will be possible. 

Consider the constraints 



C, 
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5ijd s g ks 



Using the linearized version of evolution equations 
time, 



and fll3), one can see that these constraints are constant in 



doCkij — 0, 



independently of whether or not the remaining constraints are satisfied. The constraints Ck = dk — 12dk<fi are constant 
in time as well, provided that r\ + 3x = 0. Therefore, if we give initial data with Ckij = and Ck — 0, the equations 



1 



- ^d k gij 



6 kij 

d k = 12d k 



-Is 
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1 



k{i^"9j)s + j^8ijd s gks, 



(14) 
(15) 



hold at all times, provided that »/ + 3x = 0. In this case, we can replace ekij and dk by the right-hand side of equations 
( |l4| ) and (JlJ) everywhere they appear in the linearized version of equations (0-^0 12) and forget about the evolution 



(16) 



equations for the 15 variables ekij, dk- The evolution equations for K and Aij become 
d Q K = (l + \ l\ d k T k - 8 (l + ^7 + \ A d k d k cp - e-®d k dk(e Q ), 



d k d k 



%9ii + 



% r j) + ^( 15 + 90 (d (i d s g j)s - d (i T j} ) - 2(1 + 6<r)W - e^c^) 



TF 



(17) 



In the following, we will call the system which constitutes the linearized version of equations ([7|,^|Jl2|) and the equations 
(|l6|,[l7]) system B. If initial data is chosen for system B and one takes the same initial data for system A with dk and 
ekij given by equations (^4] , [l5| ), the two systems will yield the same result upon evolution (if there are no boundaries 
or if we have time-like boundaries to which do is tangential.) 

If we compare system B with the generalization of the BSSN system by Alcubierre et al. ^ in the linearized regime, 
we see that the two agree if the KST parameters are 



c 



7 



n 



-m , x 



(18) 



Notice that since || reduces to BSSN |7|, ^ when m = 1, system B agrees with the original BSSN as well for these 
values of the parameters. Notice also that for the equivalence to be true one has to use in BSSN the same gauge 
conditions used in KST (Eq. (m) with arbitrary a). 
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For example, the choice 7 = — | eliminates the Ricci scalar appearing in the evolution equation for K, and choosing 
i] = I eliminates the term d s A s k in the evolution equation for the variable Tj, as in the traditional BSSN system 
(while in Q the parameter m is allowed to vary). Therefore, the linearized BSSN equations are dynamically equivalent 
to the linearized KST system with the choice of parameters given by (|l^) , provided that the constraints ( pT|JT5| ) are 
satisfied. In particular, any solution to the linearized BSSN equations (even one that does not satisfy the Hamiltonian 
or momentum constraint, or the constraint that corresponds to the definition of Tj) gives a solution to the linearized 
KST system with (|l^) if we define dk and e k ij according to ( [T^ , |l5| ) . 

Now an interesting question is whether or not the KST subfamily that corresponds to the linearized BSSN system 
is SH. Recalling Section II, it turns out that the characteristic speeds are plus or minus the square root of 

Ai = 2cr, A 2 = i(4m - 1), A 3 = 771, A 4 = 1. 

In particular, we see that densitizing the lapse with a > in ( |l8| ) corresponds to a SH system as long as m > 1/4. In 
particular, if a = 1/2 and m = 1, the system has only physical speeds (A = 0, 1). It is also worthwhile pointing out 
that in this (linear, constant coefficient) case well posedness can be shown for the SH families, since all is needed is 
the boundedness of the symmetrizer, which can be shown to hold. 

One might wonder whether the parameters ( [l8| ) (with a — 1/2) correspond to the ones found empirically by KST 
to have better stability properties. The answer is no, since the one-parameter subfamily ([l8]) has no intersection with 
the two-parameter Einstein-Christoffel (EC) Jl2| generalization used in the numerical simulations of KST. 



III. HYPERBOLIC FIRST ORDER ENLARGEMENTS OF THE BSSN SYSTEM 



We were not able to find a way in the non-linear case to get rid of the extra variables e k ij and dk when lower order 
terms are present in the equations. We will therefore not prove a correspondence between the KST and BSSN systems 
but, instead, a correspondence between the BSSN system and certain SH or symmetric hyperbolic formulations closely 
related to one by Friedrich and Rendall 10 . Starting with the (non-linear) BSSN system with vanishing shift (see 



below for the case of non-vanishing shift) we are able to show the following result: by introducing 18 extra variables, 
the BSSN system can be written as a system which is first order in time and space and whose principal part looks 
very similar to the one of the CT decomposition of the KST system. Drawing from our experience with the latter, 
one can then show that the new first order system has the same characteristic speeds as system B with the choice of 
parameters (|l8|). With a specific choice of parameters a and m this system has the same level of hyperbolicity as the 
KST one. 

The non-linear BSSN equations with zero shift have the form 



d (j) 
d K 

d Q T t 



1 



= -12ad k d k 4>-e- Q d k d k (e Q ) + l.o., 



-4<p 



1 r 



d k dk~ 9lJ + %r,) - 2(1 + 6a)did j( f> - e- Q didj(e Q ) 



TF 



(1.0. 



13 i 



Am , 



2(777 -\)d s A si -—d i K+(l.o.) i 



(19) 

(20) 
(21) 

(22) 
(23) 



where do = dt/N, and now "Z.o." denote terms that contain only derivatives of 4>, c/ij, Q up to first order, and no 
derivatives of the other variables. Here, the lapse N and the densitized lapse exp (Q) are related to each other as in 
equation ([[]). In terms of the BSSN variables, 

N = e 12 ^ +Q . (24) 

For generality, we have also added the momentum constraint to the evolution equation for with an arbitrary factor 
m as in Alcubierre et. al. ||. 

Our aim is to recast this system of equations into a first order system by introducing extra variables, 

d k = 12<9 fc 0, d kij = dk9ij > ( 25 ) 
and using only the constraints that are associated to the definition of these extra variables. Using the fact that 



[9o,d k 



N 



d = (12er<9 fe + d k Q) d = [ad k + d k Q) d , 
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and Eqs. ( |19| , |20] ), we obtain the following evolution equations for d k and dkij- 

d dk = -2{dk + ad k + d k Q)K, (26) 
d d Hj = ~2(d k + ad k +d k Q)A tl . (27) 

The constraints associated to the definition of d k and dkij , C k and C k ij , satisfy 

d C k = -2aKC k , (28) 
d C kij = -2aAaC h , (29) 

independently on whether or not the remaining constraints are satisfied. Since the evolution equations for C k and 
Ckij form a closed system, the equations that define d k and dkij are satisfied automatically during evolution (provided 
that they are satisfied initially and appropriate conditions are given at time-like boundaries.) 

Thus, we can enlarge system (|f9-23) by adding the evolution equations (^,^) above and by replacing 12d k <j> by 
d k and d k gij by d k ij , respectively. The resulting system is 

a o = -\k, (30) 

6 

d gij = -2Aij, (31) 
8 K = -ad k d k -e- Q d k d k (e Q )+l.o. 1 (32) 

1 - T /l \ ~\ TF 

doAij = --d k d klJ + e" 4< * d {l T j} - f - + a J d (i d j} - e- Q didj(e Q ) + (Lo.)ij , (33) 

~ 4 777 

d Tj = 2(m - l)8 s A SJ - — djK + {l,o.) 3 (34) 

d d k = -2d k K + (l.o.) k , (35) 

d d~kij = -2d k A i0 + (l.o.)kij , (36) 

where after the replacements mentioned above, no derivatives of the dynamical variables appear in "Lo.". Analyzing 
the level of hyperbolicity, it turns out that the system ( p0[ - ^6| ) has the same characteristic speeds as the KST system 
with the choice (|l8| ) and the same level of hyperbolicity. In particular, it is SH if a > and m > 1/4, and has only 
physical characteristic speeds if a — 1/2 and m = 1. 
Furthermore, it is easy to see that if the relation 

2a=i(4m-l) (37) 

holds and m > 1, the principal part of the system is symmetric with respect to the inner product associated with 

3 p-40 / Am Am 2 \ p - i( t> / \ 

(u, u) = + ,r', hl + K 2 + A^Aij + T k T k - —T k d k + —d k d k + — - (d k d k + d k ^d kij ) , 

4m — 1 2(m— 1) \ 3 9/4V / 

where u = (<j),gij, K, Aij,d k ,Tj,d k ij) T and where indices are raised with the metric gij. A point to be remarked is 
that this symmetric hyperbolic subfamily has superluminal speeds, since we have to impose the condition m > 1 and 
Eq. ([}?]) must hold. In particular, we cannot recover the original BSSN system (m = 1) but, instead, its generalization 
by Alcubierre et. al.. 

To summarize, we have shown that any version of BSSN which has the form ( |l9[|23| ) can be rewritten as a first order 
system in time and space by introducing the extra variables (p5|). Provided that initially the constraints associated 
to the definition of the new variables are satisfied, this first order system is equivalent to the BSSN system chosen. 
For a fixed densitized lapse, the system is SH if a > and m > 1/4. If m > 1 and we choose to densitize the lapse 
according to Eq. (^), the system is even symmetric hyperbolic. At this point it is worthwhile pointing out that a 
similar enlargement to a SH first order system does not work in the ADM case, even after densitizing the lapse. In 
order to get a SH system, one must add the Hamiltonian or momentum constraints to the evolution equations. 

If the shift is non-trivial, one can still obtain the above enlargement. But in this case, divergence terms of the 
conformal metric have to be substituted by T^'s in such a way that the equations have the form jl9|-p3|) with 
do = (dt — j3 k d k )/N and where the lower order terms ("Lo.") may depend on /3\ In this case, the enlargement is SH 
or symmetric hyperbolic under the same conditions as is the zero shift case provided the shift is a prescribed function. 

Frittelli and Reula (FR) have also presented well posed CT formulations of Einstein's equations [|l3| that have 
attractive properties (for example, only physical characteristic speeds), but they do not seem to reduce to the BSSN 
system in the sense described above due to the way the constraints are added in their construction. 
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IV. SUMMARY 

We have shown that at the linearized level the BSSN system with hxed densitized lapse and fixed shift is equivalent 
to a SH subfamily of the KST formulation, while at the nonlinear level it is equivalent to certain SH or symmetric 
hyperbolic system. This equivalence is analogous to the one between the wave equation written as first order in time 
and space, or first order in time but second order in space. 

The hyperbolicity of BSSN (and any other system) depends delicately on how the equations are written. Small 
substitutions of some equations in others are enough to spoil hyperbolicity. This may account for several anecdotal 
reports that rewriting part of the BSSN equations appears to improve or worsen the numerical behavior. Having 
established that this system is hyperbolic opens the possibility of rigorously showing that the resulting numerical 
schemes constructed upon discretization are stable in the sense of Lax's theorem, and therefore the resulting codes 
are convergent. This property is not available at present for WH formulations, including the ADM equations or 
even small modifications of the BSSN system considered here. Another advantage is that one could also derive 
consistent boundary conditions for the BSSN equations or show well posedness by looking at the first order hyperbolic 
enlargement and later reducing the result back to second order in space. This therefore lays mathematical grounds 
on which to understand and make use of the advantages in numerical implementations of the BSSN equations. 
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